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ABSTRACT 

We introduce the Minimum Entropy Method, a simple statistical technique for constraining the Milky Way 
gravitational potential and simultaneously testing different gravity theories directly from 6D phase-space sur- 
veys and without adopting dynamical models. We demonstrate that orbital energy distributions that are separa- 
ble (i.e. independent of position) have an associated entropy that increases under wrong assumptions about the 
gravitational potential and/or gravity theory. Of known objects, 'cold' tidal streams from low-mass progenitors 
follow orbital distributions that most nearly satisfy the condition of separability. Although the orbits of tidally 
stripped stars are perturbed by the progenitor's self-gravity, systematic variations of the energy distribution can 
be quantified in terms of the cross-entropy of individual tails, giving further sensitivity to theoretical biases in 
the host potential. The feasibility of using the Minimum Entropy Method to test a wide range of gravity theo- 
ries is illustrated by evolving restricted N-body models in a Newtonian potential and examining the changes in 
entropy introduced by Dirac, MONDian and f(R) gravity modifications. 



1. INTRODUCTION 

Studies of the Milky Way potential usually involve the con- 
struction of dynamical models under a given law of gravity 
(e.g. Newton's). Knowledge about the amount and distribu- 
tion of mass in our Galaxy is typically gained by solving the 
equations of motion in a given gravitational potential and con- 
trasting the body motions predicted by the theoretical models 
against the observed kinematics of different sorts of tracers 
(i.e. bulk gas motions; stellar radial velocities and proper mo- 
tions). 

Statistical tools provide an alternative methodology for in- 
ferring the dynamical properties of gravitational systems di- 
rectly from observational data sets and without an a priori un- 
derstanding of the motions of kinematic tracers in the system 
that is being observed. In this contribution it will be shown 
that under some special conditions it is possible to derive 
statistical inferences about the Milky Way potential without 
solving the equations that govern the motion of stars in the 
Galaxy, but by merely assuming that the orbital energy is an 
integral of motion. Gravity tests can thus be incorporated in 
the analysis in a quite straightforward manner. 

The method of statistical inference outlined here focuses 
on stellar tidal debris, i.e. stars that are tidally stripped from 
gravitationally-bound objects. These stars have a unique dy- 
namical property: although they sample very large volumes 
of the host's phase-space (position-velocity), their orbits are 
close to that of the progenitor system. The Sagittarius dwarf 
galaxy provides a dramatic example, as it is currently shed- 
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ding stars to the Milky Way halo in the form of a tidal stream 
that completely wraps around our Galaxy (e.g. Koposov et al. 
201 lb and references therein). It has been extensively shown 
that the location and kinematics of tidal tails can be used to 
put strong constraints on the host potential under the assump- 
tion that tidal tails follow single orbits (e.g. Law et al. 2005, 
2009; Penarrubia et al. 2005; Jin & Lynden-Bell 2007; Bin- 
ney 2008; Eyre & Binney 2009; Koposov et al. 2010). 

Unfortunately, most tidal debris will be too faint to be de- 
tected as coherent over-densities (e.g. tidal tails, clouds or 
shells) in photometric or spectroscopic surveys, either be- 
cause their progenitors contained a small number of stars, or 
because the stars were stripped a long time ago. Indeed, the 
dynamical evolution of tidal debris adds considerable com- 
plexity to the detection, follow-up and subsequent dynami- 
cal analysis of these systems, for tidally-stripped stars tend to 
progressively fill the allowed phase-space volume of the pro- 
genitor's orbit, becoming in the process kinematically colder 
and spatially more sparse (e.g. Helmi & White 1999). 

In contrast, the evolution of tidal debris in the space of inte- 
grals of motion is remarkably quiescent. In this space tidally- 
stripped stars are distributed tightly about the integrals of mo- 
tion of the progenitor object regardless of how far in the past 
the disruption event occurred; a property that holds even in 
potentials that vary slowly with time (Penarrubia et al. 2006). 

Some integrals of motion, like the angular momentum in 
a spherical potential, are quantities that can be measured di- 
rectly from phase-space coordinates. Others, like the orbital 
energy in a static potential, require additional hypothesis of 
the behaviour of gravity in the Milky Way. Recently, Bin- 
ney (2008) found that wrong theoretical assumptions on the 
Milky Way potential increase the scatter in the orbital ener- 
gies inferred from the locations and velocities of individual 
stars along tidal tails 5 . This statistical property enables him 

5 McMillan & Binney (2008) show that a similar principle applies in the 
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to construct an algorithm for measuring the Galactic potential 
by minimizing the r.m.s. energies derived from fitting single 
orbits to tidal tails. However, this algorithm relies on the as- 
sumption that all stars on the tails have the same orbital energy 
in the true gravitational potential. 

Here it is demonstrated that Binney's (2008) results can be 
extended to tidal debris with arbitrary energy distributions as 
long as these distributions are separable (i.e. invariant) in 
space. This assumption may be accurate for stellar streams 
whose progenitors have low dynamical masses, the so-called 
'cold' tidal streams (Penarrubia et al. 2006; Eyre & Bin- 
ney 2009). Also, it will be shown that these systems can be 
used to measure the Milky Way gravitational potential as well 
as to test different gravity theories simultaneously, and with- 
out having to integrate the equations of motion. Indeed, our 
method provides a simple tool to constrain the Milky Way 
gravity through a direct inspection of the phase-space coordi- 
nates of cold tidal streams, thereby avoiding the daunting task 
of solving for the dynamics of stars in host galaxies with a 
generic mass distribution under a suite of gravity theories. 

Dropping the necessity to fit orbits to the locations and ve- 
locities of streams brings a further advantage, namely the fact 
that our method can be equally applied to phase-space sub- 
structures that can be identified as tidal tails, as well as to sub- 
structures that have 'dissolved' in the Galactic halo, but which 
may be detected as substructures in the space of integrals of 
motion, as long as they originate from low-mass progenitors. 

Before introducing our method we wish to remark that 
the technique outlined in this contribution requires full 6D 
phase-space information. At present most stars with mea- 
sured position-velocity vectors are located in the solar neigh- 
bourhood, which severely limits its applicability. However, 
this unfortunate situation will likely improve in the near fu- 
ture thanks to the advent of Gaia, an European mission that 
will provide astrometric data for more than one billion stars 
with unprecedented accuracy (Perry man et al. 2001). Gaia 
data in combination with complementary spectroscopic sur- 
veys (e.g. Gaia-ESO) is expected to reveal a large number of 
substructures associated with past accretion events (Brown et 
al. 2005; Mateu et al. 2011). 

The paper arrangement is as follows: we develop the core 
idea of our method in §2. Tests are shown in §3. We go 
through a brief discussion about applications and limitations 
in §4 and §5. 

2. THE ENTROPY OF TIDAL DEBRIS 

In statistical mechanics, entropy measures the degree to 
which the probability of the system is spread out over different 
possible states. The concept of entropy has been widely used 
in Physics (e.g. Wehrl 1978; Jaynes 1957) to determine the 
behaviour of macroscopic systems in equilibrium, or close to 
equilibrium. For example, galaxies undergoing phase-space 
mixing tend to evolve toward an energy configuration that 
maximizes entropy 6 . 

action-angle space, where wrong assumptions on the local potential alter the 
characteristic periodicity (or "beats") of the patterns shown by tidal debris as 
they cross the solar neighbourhood. 

6 More specifically, Tremaine, Henon & Lynden-Bell (1986) show that 
the equilibrium configuration is that which maximizes all H-functions H = 



In this work we are interested in stellar systems that follow 
a narrow energy distribution and have, by construction, low 
entropy. Stars that are tidally stripped from gravitationally- 
bound, low-mass objects and orbit about the host galaxy 
potential provide an example of systems with low entropy. 
Given that biases in the theoretical modelling of our Galaxy 
tend to increase the range of energies sampled by tidal debris 
(Binney 2008), it appears therefore natural to use entropy as 
statistically-meaningful quantity to identify those biases. In 
addition, entropy has special properties that will help us be- 
low to tackle difficult aspects of the study. For example, its 
additivity, i.e the fact that the combined entropy of two in- 
dependent systems equals the sum of the individual entropies 
(see §II.E of Wehrl 1978), helps to combine information on 
the Milky Way potential gathered from individual debris sys- 
tems (§2.2), as well as to analyze the effects of the progeni- 
tor's self-gravity (§4). 

It is illustrative to re-interpret Binney (2008)'s results in 
terms of entropy. Consider a suite of TV* stars orbiting in a 
potential $(r) with phase-space coordinates {r^v,}^ and an 
orbital energy per unit mass £, = E for all i, where 

E, = |+$(r ( ). (1) 

Given that all stars have the same energy in the true poten- 
tial, the energy distribution follows a Dirac's delta function 
per construction. Binney (2008) shows that wrong assump- 
tions about <l>(r) must necessarily yield a scattered sample of 
energies unless all the stars are located at the same position. 
The larger the departure from the true underlying potential, 
the higher the rms variation in energy. Using entropy one 
would conclude that, since by definition a delta function has 
an entropy of minus infinity, the entropy of the biased energy 
distribution can only increase owing to poor choices in the 
potential model. Clearly, the stronger the bias the larger the 
scatter, and thus the higher the measured entropy must be. 

In §2.1 we extend this principle to stellar systems whose 
energy distribution f(E) is differentiable and separable in en- 
ergy and space 7 . Mathematically this condition implies that 
the probability that a star having an energy £ at a location r 
can be written as p(E,r) = f(E)g(r), where g(r) is the proba- 
bility that a star has position r. Both probability functions are 
normalized so that / f(E)dE = J g(r)d 3 r = 1. It is convenient 
to define the relative potential as ^ = -$ + $oo, and the rela- 
tive energy as e = -E + Q^. Here is an arbitrary constant 
that sets the boundary conditions for the energy distribution 
/(£)■ 

The fact that the Milky Way potential is unknown intro- 
duces a model-dependent bias ((5$) in the orbital energy de- 
rived from phase-space measurements 

e(r) = £ (r) + (5$(r). (2) 

In what follows we use tildes to denote measured quantities, 
whereas symbols without tildes denote true (i.e. unbiased) 

- / C(F)drdv, where F(r, v) is the coarse-grained distribution function, and 
C(F) is any convex function with C(0) = 0. Entropy corresponds to the case 
C(F) = F\nF 

1 Note that a similar analysis can be developed in the space of actions. 
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quantities. 



2.1. Single component 



For systems with separable energy distributions, the bi- 
ased energy distribution relates to the true one as p(e,r) = 
p[s-5$(r),r] = /[e-£<&(r)]g(r). Let us expand the measured 
energy distribution /(e) at order C[((5$) 2 ] 



f(e) 



I 



<S$(r)]g(r)d 3 r f 



(3) 



l-^(r)^ + ^ 2(r)/ " (£) 



g(r)d 3 r = 



i- m m + {6 * 2)f " i£) 



Here brackets denote volume-averaged quantities, i.e. (x) = 
J xg{r)d?r, whereas /' = d//de and /" = d 2 //de 2 . 

The entropy of the measured energy distribution is by defi- 
nition 



de/(e)ln[/(e)]. 



(4) 



After approximating ln(l +x) w x-x 2 /2 forx <C 1 and some 
algebra, the combination of eq. (3) and (4) at 0[(5$) 2 ] results 



~ff+<5$> J de/(e)[l+ln/(eXP) 



(6<S>) 2 



de/(e) 



/'(e) 
/(£). 



(<5$ 2 



de/"(e)[l+ln/(e)]. 



The integrals including the term (1+ln/) can be eas- 
ily solved by parts. Adopting boundary conditions so that 
lim/ln/= lim /In/ = lim /'In/ = lim f'lnf = 0one 

can readily show that the first term is zero 

j d £ /'(l+ln/)=(/ln/)*-=0, 



whereas 



7* 



-,2 



Hence eq. (5) becomes 



/ 



de/(e) 



fie) 

m 



where the quantity a\ = ((5$ 2 ) - (<J$) 2 is the dispersion in 
the energy distribution that arises from the bias in our Galaxy 
potential (eq. 2) over the volume occupied by the tidal de- 
bris sample, and a~ 2 = J de/(e)[/'(e)//(e)] 2 is related to the 
internal dispersion of the energy distribution in the unbiased 
Galaxy potential. For example, for a Gaussian distribution 
/(e) = 1 / V27rcr 2 exp[-£ 2 /(2(T 2 )] the correspondence is direct, 
as a~ 2 = J de/(e)[/'(e)//(e)] 2 = a~ 2 . 



Eq. (6) shows a few interesting points. First, since both 
quantities o\ and a 2 are positive, we find that any bias intro- 
duced in our Galaxy model yields an increase in the entropy 
of the energy distribution. Second, the choice of $oo = const 
does not alter the value of the entropy, as adding a constant 
bias to the potential leaves the dispersion of orbital energies 
unchanged. Third, because AH ~ l/2(<r^/a E ) 2 we find that 
the change in entropy is fairly sensitive to biases in the poten- 
tial averaged over the region probed by the tracer population. 
This is a remarkable property of the entropy, suggesting that 
minimization of entropy for stellar systems with separable en- 
ergy distributions provides a powerful method to measure the 
Milky Way potential. Cold tidal debris from low-mass pro- 
genitors may be an example of such systems because they oc- 
cupy a reduced volume in the integral-of-motion space, but 
they sample large volumes of phase-space (see §4). 

2.2. Multiple components 

In practice it may prove challenging to isolate the energy 
distribution of a single population of stellar debris from the 
Milky Way background, or from a population of distinct struc- 
tures overlapping in a given energy range. It is thus instruc- 
tive to repeat the above calculations for an energy distribution 
of the form f(e) = a/i(e) + (l -a)/ 2 (e), where < a < 1. 
The measured distribution is thus /(e,r) = a/i(e)gi(r) + (l- 
a)f 2 (e)g 2 (r). 

After some algebra it is straightforward to show that a bias 
in the galaxy potential of the form given in eq. (2) leads to a 
change in the entropy 

Aff ~a<<y*i) J de/'tl+ln/J-a^ J de/ftl+ln/] 



-a 



. (1 . a) M)|^ [1+hfl . (1 . a) ^ J def 



-a(l-a)(S^)(S^ 2 ) / d£ ^) 



Note that now the terms that contain the linear variations 
of the potential, (<5$), in general do not vanish. Hence the 
entropy of the function /(e) may increase or decrease de- 
pending on the potential bias and our choice of $oo. These 
terms only vanish if (i) the structures defined by /(e) and 
/2(e) follow the same spatial distribution, i.e. gi(r) = #2(1") 
so that (<5$i) = (<5<i>2) and (S$ 2 ) = {5<&\); or (ii) if the en- 
ergy distributions of the substructures do not overlap, so 
that / def [ 1 + In /] = / A£i de/' [ 1 + In /•] ; and / de/(/' / f) 2 = 

Iac d £ Mf! I ' fd 2 , where Ac, is the range of energies occupied 
by the i'th substructure. In this case the crossed terms van- 
ish, J de(/'/j)// = 0, and by analogy with §2.1 we find the 
entropy variation of N s non-overlapping structures is AH w 
1 /2Y^1\ Oii(a^j/a £ j) 2 , where Y^,f=i a i = 1- This is indeed an 
interesting result, for it allows us to combine in a simple way 
the constraints on the Milky Way potential derived from dis- 
tinct substructures under the only condition that their energy 
distributions do not overlap with each other. 
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FIG. 1. — Left column: Spatial distribution of particles that follow Gaussian orbital energy distribution with a dispersion <j £ = 10 -3 for orbits with three 
different orbital pericentres (rows). The orbital apocentre of these models is 5do. Middle column: Effects on the orbital energy distribution of introducing a 
bias in the value of $0. The un-biased distribution is Gaussian (solid lines), whereas dotted and dashed lines denote biases of +1% and —1% in the value of <E>o. 
respectively. Note that the changes in the measured energy distribution become more prominent as the orbital eccentricity increases. Right column: Effects of a 
bias in the value of do. As expected, the changes in the energy distribution are small in models where all particles are located at r S> do. 



Finally, eq. (7) shows that, as in many other problems in 
Physics, it will be necessary to clean out the sample of stel- 
lar tidal debris from Milky Way background objects, as that 
is a case of two overlapping distributions by definition. To 
this end it may help to consider additional integrals of motion 
(e.g. the vertical component of the angular momentum if the 
potential is axi-symmetric), and/or metal compositions. Also, 
at some point it may be worth extending our analysis to the 
space of actions (e.g. McMillan & Binney 2008), where the 
background subtraction may be more straightforward. 



3. TESTS 

In this Section we devise a number of numerical experi- 
ments that aim to test the analytical results enclosed in §2. 
For simplicity we adopt an unbiased energy distribution that 
is Gaussian, /(e) = l/^/27raf exp[— (e- £ or b) 2 /(2o^)], where 
e or b is the orbital energy of the disrupted system and a e its 
dispersion. Recall, however, that the results obtained in §2 
apply to any energy distribution that is differentiable in the 
range [0, $oo] and separable in energy and space. 
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FIG. 2. — Entropy of the orbital energy distribution as a function of bias in 
the value of $o (left column; Sdo = 0) and do (right column; 8&o = 0) for 
three different orbits. The entropy of the unbiased (5$o = Sdo = 0) Gaussian 
distribution is H = 1 /2[ln(27rcr2) + 1] ~ -5.49 for a e = I0r 3 , and H ~ -7.79 



for <j e 



10" 



Note the strong sensitivity of the entropy to biases in the 



potential parameters as <5<3?o and Sdo approach zero (i.e. unbiased potential). 
3.1. Set-up 

Again for simplicity we adopt a host galaxy potential that 
is spherical and does not vary with time 



$(r) = $ ln(d +'- 2 ). 



(8) 



Subsequently, we construct suites of 10 4 particles ini- 
tially placed at an initial radius rn. The initial velocities 
v = \/2[— $(ro) — e] are chosen so that the energy distribu- 
tion is Gaussian. The mean of the energy distribution is 
£ rb = -Vo/2-$(r ), where v = £v c = £y/rd$/dr, and v c is 
the circular velocity of the host. The parameter £ < 1, so that 
all our orbits are initially at apocentre. 

The initial azimuthal location of the particles is random. 
The integration time falls between 1 and 10 t a , where t a = 
ro/vo. The equations of motion are solved using a leap-frog 
scheme with a time-step that is chosen to provide an energy 
conservation AejS 10 -1 a e . 

Finally, throughout this Section we adopt the N-body units 
G= <f> = d = l. 

3.2. Models 

For illustrative purposes all the orbits considered here have 
an apocentre ro = 5do = 5 (we note, however, that a different 
choice of ro would not change the qualitative results shown 
below). By varying the parameter £ one can change the or- 
bital eccentricity of our models in a simple way. Values of 
£ = and 1 yield radial and circular orbits, respectively. To 
illustrate the method outlined in §2 we consider models with 
£ = 0.1,0.6 and 0.9, which respectively yield orbital pericen- 
tres 0.27 (orbit C), 2.14 (orbit B) and 4.11 (orbit A). Their 
mean orbital energy is e or b ~ -4. 1 1 ,-3.60 and -3.27. 



We choose energy dispersion values for our models that 
are are representative of those expected for the tidal debris 
of dwarf spheroidal galaxies. The dispersion of the orbital 
energy distribution relative to the Galaxy potential roughly 
scales as a e ~ Cv/ V max< where a v is the original velocity 
dispersion of the disrupted system and v max is the peak ve- 
locity of the Milky Way rotation curve. Adopting v max = 
220kms _1 we expect a e ~ 10 -4 — 10~ 3 for dwarf spheroidal 
galaxies, which show velocity dispersions in the range ~ 3 
and lOkms -1 , respectively (e.g. Mateo 1998; Simon & Geha 
2007; Koposov et al. 2011a). Although low-mass globular 
clusters can have velocity dispersions as low as a v ~ lkms , 
which correspond to a e ~ 10~ 5 , here we only consider inter- 
mediate values for the energy dispersion, a e = 10 -4 and 10~ 3 . 

3.3. Biases arising from the potential parameters 

The spatial distribution of test particles for the three orbits 
considered here is plotted in the left column of Fig. 1. As 
shown by the solid lines in the middle and right columns, 
these particles follow a Gaussian energy distribution in the 
unbiased (6&o = Sdo = 0) potential. The middle column illus- 
trates the effect of introducing a small bias in the value of <£>o 
on the shape of the energy distribution. Note that biases in 
the potential change both the mean orbital energy (e) and the 
energy dispersion of the particles. However, by plotting the 
energy distribution as a function of (e— (e)) we can only ap- 
preciate variations in the latter. Red-dotted and blue-dashed 
lines correspond to values of <E>o = 1.01$o an d $o = 0.99$o, 
respectively. As one would expect, the shape of the energy 
distribution changes more strongly as the orbital eccentricity 
increases, and thus the range of orbital radii, increase. The 
right column shows the effects of changing do. Here it is inter- 
esting to observe that the energy distribution barely changes if 
all particles move on orbits with r^ do (orbit A), highlight- 
ing the fact that constraints on a given potential parameter can 
only be derived from orbits that are sensitive to variations in 
that particular parameter. Indeed, as more particles are lo- 
cated at r ~ do (e.g. orbit C) changes in f(e) become again 
visible. 

The entropy associated to the orbital energy distribution of 
these models (eq. 4) is calculated using a simple trapezoidal 
rule. The results are plotted as a function of biases in $q and 
do in the left and right columns of Fig. 2, respectively. Rows 
correspond to models with an energy dispersion a e = 10~ 4 
and 10~ 3 . This Figure illustrates two interesting points. First, 
the entropy is indeed minimum for the unbiased potential, as 
analytically predicted in §2. Biases in any of the potential 
parameters translate into an increase of entropy that is more 
marked for models with intrinsically narrow energy distribu- 
tions. Second, the increase in entropy becomes more pro- 
nounced as (5<&o and Sdo approach zero. This is an important 
result, as it suggests that entropy minimization may provide 
a powerful method to measure the gravitational potential of 
galaxies with high accuracy. 

Fig. 3 shows the relative variation of entropy AH as a func- 
tion of <7$/c7 £ , where o$ = yj (5$ 2 ) - (5$) 2 is the disper- 
sion in the potential introduced by a model bias. For small 
potential biases (o$ <C er e ) the entropy increase scales as 




FIG. 3 . — Variation of entropy as a function of the dispersion in the potential 
values (cj.) introduced by biases in <J>o (similar curves are obtained from 
biases in do). The dashed line shows the theoretical expectation from eq. (6). 
Note that for small biases in the potential the entropy of the orbital energy 
distribution goes as AH ~ 1/2(<j$ /<t e ) 2 . 

AH = 1 /2(cr$ / a e ) 2 independently of orbital eccentricity and 
energy dispersion, as expected from eq. (6). For cr$ ^S> o E 
the Taylor expansion used in §2 does not hold, and deviations 
between the analytical expectation and the numerical results 
start to arise. These results provide the basis of our method 
for measuring the Milky Way potential, as they show that the 
true Milky Way potential is that which minimizes the entropy 
of cold tidal debris with separable energy distributions. 

3.4. Biases arising from the potential form 

Biases may also arise from a wrong parametrization of the 
Milky Way potential. This is illustrated in Fig. 4 by means of 
a simple experiment. Let us consider the following potential 



FIG. 4. — Effects of a bias in the functional form of the potential model. 



$(r) = 2$ 



y + y - + . 
y 3 



(N-l)/2 y2k+l 



k=0 



2k+l 



(9) 



where y = (r/do) 2 /[2 + (r/do) 2 ]. One can easily show that 
limjv->oo & = 3>, hence we expect a minimum in the entropy 
in that limit. This is clearly seen in Fig. 4, which shows an 
entropy that decreases toward the entropy of the Gaussian dis- 
tribution (dotted lines) as the number of terms in the potential 
series increases. Thus, this exercise highlights the possibility 
of using the entropy of stellar tidal debris not only to find the 
best parameters of a Galaxy potential, but also to distinguish 
between different parametrizations. 

3.5. Biases arising from the adopted law of gravity 

Let us now consider other gravity theories that allow for 
the definition of an integral of motion equivalent to the New- 
tonian orbital energy. From §2 it follows that choosing the 
wrong gravity or cosmology shall introduce a bias in the or- 



Here we expand the host potential model as <J>(r) = 2<J?o 2~Zl-o 



CV-D/2 ,,2/Sr+l 



/(2k + 



1)1 + lnrfj), with y = (r/do) 2 /[2 + (r/do) 2 ]; hence the true logarithmic po- 
tential of eq. (8) is recovered in the limit N — > oo. For ease of reference we 
plot dotted lines to mark the entropy of a Gaussian distribution. As expected, 
the entropy asymptotically converges to its minimum as the number of terms 
in the potential expansion increases. 

bital energy values, which is bound to increase the entropy 
associated to the energy distribution of tidal debris. 

3.5.1. Time-variability of G 

Let us examine first a cosmology where gravity is Newto- 
nian, but G evolves with time. For example under Dirac's 
large-number hypothesis (Dirac 1938) 8 

„2 



Gm p m e 



10 



-39 



(10) 



where t is the time since the big bang, m p and m e are the 
proton and electron masses, and e is the electron charge. In a 
Universe where the properties of elementary particles remain 
constant, G oc l/t. 

By inspection of the equations of motions in a Dirac cos- 
mology, Lynden-Bell (1982) found that the orbital energy per 
unit mass of a particle moving in a potential $ associated to a 
mass distribution p is a constant of motion if written as 

2 



■.Hit 1 



drY 
dt I 



+ ^-$(r)- 



dr 

dt 



// 2 r 2 ;(ll) 



where V 2 $ = 4-irGop, Hq is the Hubble constant and G = 
Go /(Hot). Here the sub-index "D" denotes quantities calcu- 
lated in a Universe where G oc l/t. 

Thus, by analogy with eq. (2) the energy bias introduced 
by the cosmological model at the present time, t = Hq , is 
S<Z> D = ±[-H (dr/dt ■ r/t)+l /2// 2 r 2 ], where the minus and 
the plus symbols correspond to Universes with a constant and 

8 See Uzan (2003) for a compilation of experimental bounds on the varia- 
tion of G with time. 
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an evolving G, respectively. Two points are noteworthy: first, 
S$ D only vanishes in the case of centrally-located (r = 0) or- 
bits; in the case of circular orbits, for which dr/dt • r = 0, it 
reduces to 1 /2HqT 2 . And second, the energy bias introduced 
by our choice of gravitational potential (see §3.3 and §3.4) is 
independent of that associated to our choice of cosmology. In 
practical terms this implies that the Milky Way potential and 
a possible time-variability of G can be constrained simultane- 
ously by minimizing the entropy of stellar tidal debris. 

To illustrate this point let us re-calculate the entropy of the 
models introduced in §3.2 using eq. (11). The presence of 
Ho introduces a physical scale in the solution, which forces 
us to scale the test-particle models to physical units. Given 
that our method will be mostly applied to Milky Way objects, 
we choose <£>o = 1 /2(220kms -1 ) 2 , do = 12 kpc and t = Hq = 
14Gyr. A comparison between the Newtonian and Dirac en- 
tropy variation as a function of biases in $o and do are shown 
in Fig. 5. For this particular example the test-particle models 
are calculated in a Newtonian potential, so the case G = const, 
(solid lines) yields AH = for J$ = M) = 0. Dotted lines 
show the effects of assuming G oc I ft. As expected, the en- 
tropy is higher than in the Newtonian case independently of 
the assumed values for <£>o and do. Note however the effects 
of a time-variable G are mostly visible at <5<I>o rj Sdo « 0. 

3.5.2. MONO 

The results obtained in the above Section also apply to grav- 
ity theories that are not Newtonian. The modified Newtonian 
gravity QMOND (Milgrom 2010) is an interesting example. 
In this theory the gravitational force per unit mass can be sim- 
ply written as 



gM = gNV(f) = g A 



1 




(12) 



where a§ « 1.2 x 10~ 10 m/s 2 is Milgrom's constant, and g N = 
-GM{< r)/r 2 is modulus of the Newtonian specific force. In 
the deep-MOND regime, g^ -C flo> we have gM s» ^JaogN', 
wheres if g^ 3> ao the Newtonian solution is recovered. In 
systems with spherical symmetry the corresponding gravita- 
tional potential is 



(13) 



whereas the mass profile associated to the logarithmic poten- 
tial of eq. (8) can be written as 

-3 



M(< r) = 



2$n 



G r 2 +d 2 ' 



(14) 



Clearly, the Newtonian acceleration reaches a maximum at 
r = do. For the fiducial parameters $o = l/2(220kms -1 ) 2 
and do = 12 kpc, we have that gN/ao( r = do) = ^o/(^o«o) — 
0.544; hence the minimum Mondian-to-Newtonian ratio is 
v(r = do)~ 1.944. 

The differences between $m and $ai are thus much stronger 
than those introduced by a time-varying G (§3.5.1). A com- 
parison between the dashed and solid lines in Fig. 5 shows 
that AH is dominated by the modification of the Newtonian 




FIG. 5. — Increase in entropy due to biases in the potential parameters 
$o> <kt> as we U as m the gravity theory, for a model with an unbiased en- 
ergy dispersion cr £ = ICT 4 ^. Orbits are integrated in a Newtonian potential 
(solid lines). As expected, the entropy finds a minimum (AH = 0) when the 
biases in the potential parameters vanish, i.e. 5<E>o = 8do = 0. To allow a 
comparison between different gravity theories we scale these parameters to 
<£>o = (220kms~') 2 /2 and cIq = 12 kpc, so that orbital energies are measured 
in units of (kms -1 ) 2 . Dotted and dashed lines show the effects of assum- 
ing Dirac's cosmology (i.e. G oc 1/f) and a MONDian gravity, respectively. 
Long-dashed lines show the same calculation in a f(R) gravity theory with 
/3 = 0.82 (n ^ 3.5) and r c = d (see §3.5.3). 

gravity given by eq. (12), rather than by the biases in the po- 
tential parameters. It is interesting to notice that the minimum 
entropy occurs for <5$o — -0.8 < 0, as one would expect given 
that the MOND potential is stronger than the Newtonian one 
at all radii. 

3.5.3. f(R) gravities 

Another interesting case of gravity models that aim to mod- 
ify Einstein's General Relativity at large scales can be found 
in f(R) theories. In these models the action can be written as 



A = 



d 4 x 



£[/(*) + £,„]; 



(15) 



here f(R) is a generic function of the Ricci scalar curvature 
that is usually assumed to follow a power-law function, i.e. 
f(R) = foR"', gfw is the metric, and C m is the standard matter 
Lagrangian. Einstein's General Relativity with contribution 
of a cosmological constant is recovered for n = 1 and f(R) = 
R + 2A. 

Capozziello et al. (2007) explore modifications of the New- 
tonian gravity where the potential of a spherical shell of mass 
dm can be written as 



d$: 



Gdm 



1 + 



(16) 



where 



12« 2 - In - 1 - V36n 4 + 12n 3 - 83« 2 + 50n + 1 
6« 2 -4n + 2 
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so that the relativistic case n = 1 corresponds to /? = 0. Cases 
where 1 - (3 > yield an increase in the gravitational force on 
scales r^r c . 

The potential of an extended distribution of mass with a 
density profile p(r) can be thus written as $r = 1 /2(<f>iv + $ c ), 
where & N is the Newtonian potential and 



$ c (r) = -4ttG 



1 / dr'p(r'y 
r Jo 

i: 



(17) 



dr'p(r')r' 



In contrast to (3, the quantity r c is not universal and has to 
be fitted on individual galaxies. Interestingly, Capozziello et 
al. (2007) find that n = 3.5 (/3 ~ 0.82) provides a good fit to 
the rotation curves of several low-surface-brightness galaxies 
as well as to Hubble diagram derived from Type la supernova 
with no dark matter. 

In order to illustrate how entropy may help to constrain the 
value of (3 from Milky Way phase-space surveys we plot with 
long-dashed lines in Fig. 5 the variation of entropy as a func- 
tion of biases in $o and do for j3 = 0.82 and r c = do- Here we 
consider test-particle models evolved in a Newtonian gravity 
with an unbiased energy dispersion a £ = 10 _4< E»o- Models are 
scaled to the following physical units, $o = (220kms _1 ) 2 /2 
and do = 12 kpc, and energies are measured in units of 
(kms -1 ) 2 . As with the MOND theory, the fact that the mod- 
ified gravity is stronger than the Newtonian one at all radii 
leads to an entropy that is minimum at <5$ — ~0-4 < 0, and 
to a scale do that is poorly constrained. 

As a final remark we note that similar restricted N-body ex- 
periments can be carried in Dirac, MONDian and f(R) gravi- 
ties, where for a given mass distribution the minimum entropy 
associated to the true gravity theory must be lower than in the 
Newtonian case independently of the potential parametriza- 
tion. 

4. EFFECTS OF SELF-GRAVITY 

Tidal streams formed in realistic galaxy potentials are 
poorly represented by single orbits (e.g. Penarrubia et al. 
2006; Eyre & Binney 2011). Numerical experiments show 
that unbound particles escape through the Lagrange LI and 
L2 points of the disrupting object (e.g. Renaud et al. 2011), 
forming leading and trailing tails with orbital energies that 
are respectively higher and lower than that of the progenitor 
system. Given that tidal tails have overlapping energy distri- 
butions (Penarrubia et al. 2006) and occupy different volumes 
in phase-space, it follows from §2.2 that the progenitor's self- 
gravity must introduce a bias in our method. 

To determine the extent to which the leading and trailing 
tails have different energy distributions it is useful to measure 
the Kullback-Leibler (1951) divergence of the individual tails, 
which is defined as 



fi{e) In 



de = -Hi+H c 



(18) 



./(£). 

where the sub-index i = l,t denotes leading and trailing 
tail distributions, respectively. The right-hand term H c i = 
-J fi(s)ln f(e)de is called crossed entropy. 



In the case of a separable energy distribution one has that 
fi(e) = f(e), so that the crossed entropy of the tails is equal 
to the overall entropy, i.e H C j = H; hence, the KL-divergence 
is Di = 0. The Kullback-Leibler (or KL) divergence provides 
therefore a measure of the separability of the orbital energy 
distribution. 

Using the notation of §2.2 we have that / = afi + (1 — oi)f t , 
so that the overall entropy can be written as 



H = -J f(e)\nf(e)de = (19) 

-a J f,(s)lnf(s)ds-(l-a) J f t (e)\nf(e)de = 
otHx + (1 - a)H, + aDi + ( 1 - a)D, = (H) h ,+ (D) u ; 

where (#)/,< and (D)ij are the combined entropy and KL- 
divergence of the tidal tails. 

To illustrate how the KL-divergence varies as a function of 
biases in the calculus of orbital energy we run self-consistent 
N-body simulations of tidally-disrupting clusters moving on 
the orbit C in the potential defined by eq. (8). For ease of ref- 
erence the potential parameters are set to $o = (220kms -1 ) 2 /2 
and do = 12kpc. The cluster models follow a cored Dehnen 
(1993) profile, p c = (3M/4n)a/(a + r) 4 , with masses M = 
10 5 , 10 6 and 1O 7 M and radii a = 0.29,0.62 and 1.33 kpc, re- 
spectively, which are therefore comparable to the dynamical 
masses and sizes of dwarf spheroidal galaxies (e.g. Walker et 
al. 2009). The tidal evolution of the clusters is followed us- 
ing SUPERBOX (see Fellhauer et al. 2000 for details) for 7.35 
Gyr. The numerical resolution and time-steps are chosen so 
that over that time interval the energy is conserved at a level 
As^0.01GM/a for models evolved in isolation. 

The left panel of Fig. 6 shows the spatial distribution of 
debris projected on the progenitor's orbital plane. Here par- 
ticles are coloured-coded according to their orbital energies: 
red (leading) and blue (trailing) particles denote energies that 
are lower and higher, respectively, than that of the progenitor 
system. In all models the number of particles in the leading 
and trailing tails is practically the same, so a« 1/2. Note 
that for a fixed orbit and integration time it takes many more 
orbits for cold streams to spread a substantial fraction of 2ir 
in orbital phase. 

The middle and right panels show the variation in entropy 
and KL-divergence as a function of biases in $o and do- A 
few interesting points are worth noticing. First, in contrast to 
the tails entropy (Hi) the KL-divergence (D,) finds a maximum 
for the unbiased potential (5$ = 0. For energies measured in 
units of (kms -1 ) 2 we find that the maximum combined KL- 
divergence is (D)ij ~ 0.65, a value that barely depends on the 
progenitor mass. In contrast the combined entropy of the tails 
in the unbiased potential is (H)i, t ~ 7.17,7.95 and 8.72 for 
M= 10 5 , 10 6 and 10 7 M Q , respectively. As expected, the ratio 
between the combined KL-divergence and the entropy of the 
tails ((D)ij/(H)i, t ) is small, which indicates that tidal debris 
of low-mass objects follow nearly-separable energy distribu- 
tions. However, it is worth noting that although this quan- 
tity becomes smaller as the progenitor mass decreases, it only 
vanishes in the limit M — > 0. 
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FIG. 6. — Left panel: Spatial distribution of tidal debris. Particles are coloured-coded depending on whether they belong to the leading (in red) and trailing (in 
blue) tails. Dotted lines mark the radius r = cIq. Middle and right panels: Variation of the debris entropy (H, solid lines) as a function of biases in the values 
of 4>o and do, respectively. Energies are measured in units of (kms~') 2 . The entropy (//,) and the KL-divergence (Df) of the individual tails are shown with 
dotted and dashed lines, respectively. As expected from eq. (6) the entropy of the individual tails increases under the presence of biases in the calculus of the host 
potential. Note that the opposite trend is visible for the variation of D;, which is maximum for the unbiased <5<J>o = <Wo = potential. 



Note also that if each of the tails followed a separable en- 
ergy distribution, then from §2.1 one should expect AH, ~ 
((7$ i/tT E I ) 2 in the limit <5<!>/<P <C 1. However, this is clearly 
not the case. The reason can be found in the time-variation 
of the progenitor's self-gravity. As mass stripping progresses 
the location of the Lagrange points LI and L2 move closer to 
the centre of the disrupting system, which introduces a depen- 
dence between the energy of the particles and the time when 
they become unbound. This inevitable effect is responsible 
for the asymmetric (and mirrored) dependence of Hi and H, 
as a function of <5<Po and Sdo, which is not visible in Fig. 2. 

Finally, the fact that (!>)/., ^ as well as a decreasing func- 



tion of the energy bias has two undesirable effects. First, H is 
less sensitive to the presence of biases in the energy calcula- 
tion than (H)i ,. Second, the condition for minimum entropy 
H' = {H)' lt + (£>); , = 0, where /' = d//d<S$, is now met for 
(H)' lt = (D)' lt = 0, which yields a local maximum at <55> ss 0; 
and for |(ff)J,| = — |(D)J,|, which yields two local minima 
around <5$ « 0. The progenitor's self-gravity therefore in- 
troduces a bias in the Minimum Entropy Method. For the 
N-body satellite models explored here, which have dynamical 
masses M^S 1O 7 M and follow eccentric orbits about a Milky 
Way-like potential, the bias is fairly small, |(5$| j50.03$o- 
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However, notice that the bias does not only depend on the 
progenitor's self-gravity. The fact that AH t ~ (cr<s>.i/<Je,d 2 at 
5$ « implies that the steepness of the function (H)' l t , and 
thus the location of the two minima, will also depend on the 
orbit of the progenitor about the host potential that is being 
inspected. 

These results suggest that the effects introduced by the pro- 
genitor's self-gravity can be lessened in two ways. First, one 
may attempt to minimize the ratio (D)ij/\(H)ij\ by applying 
the Minimum Entropy Method to tidal debris that originate 
from low-mass systems. In addition, one may derive the prob- 
ability that a given star belongs either to the leading or trailing 
tail by inspection of the distribution of debris in the integral- 
of-motion space. As shown in Fig. 6, knowing the member- 
ship probability allows a derivation of the crossed-entropy of 
the individual tails, giving further sensitivity to biases in the 
host potential. 



5. SUMMARY AND DISCUSSION 

In this contribution we have introduced the Minimum En- 
tropy Method. This statistical technique is devised to con- 
strain the Milky Way gravitational potential and test differ- 
ent gravity theories directly from stellar 6D phase-space cata- 
logues and without adopting dynamical models. Our method 
rests upon two fundamental assumptions: (1) the gravity the- 
ory under study allows for the existence of an orbital energy 
that is an integral of motion; and (2) cold structures in phase- 
space resulting from the tidal disruption of gravitationally- 
bound systems have an orbital energy distribution that is sep- 
arable in energy and space. Then it is shown in §2 that any 
bias in the calculus of the orbital energy (due to the adoption 
of an incorrect potential and/or gravity theory) translates into 
an increase of the debris entropy. Examining what type of 
gravity theories obey the first condition goes beyond our cur- 
rent goals, but it is worth considering here in what cases the 
second condition may not be met. 

In §4 we show that the progenitor's self-gravity induces the 
formation of leading and trailing tidal tails, which occupy dif- 
ferent phase-space volumes and follow distinct energy distri- 
butions. Because the separability condition is broken, self- 
gravity introduces a bias in the location of the minimum en- 
tropy, thus limiting the application of the method to tidal de- 
bris that originate from low-mass systems. However, the ef- 
fects of self-gravity can be lessened by deriving the probabil- 
ity that a star in the debris sample belongs either to the lead- 
ing or trailing tail. With this information at hand it is possible 
to quantify systematic variations of the energy distribution in 
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terms of the cross-entropy of individual tails and subtract the 
contribution of the non-separable term ( (the so-called 
KL-divergence) from the overall entropy H. The resulting 
quantity is the combined entropy of the leading and trailing 
tails (H)ij = H- (D)ij, which is more sensitive to biases in 
the host gravity than the overall entropy. 

Perturbations on the orbit of a system undergoing tidal dis- 
ruption may also break the separability of the energy distri- 
bution. For example, a drag force (e.g. dynamical friction) 
acting upon a tidally-disrupting system introduces a depen- 
dence between the time when stars are lost to tides and their 
mean orbital energy. Similarly, interactions with bound sub- 
structures lingering in the Milky Way halo (e.g. molecular 
clouds, stellar clusters) are bound to introduce similar biases 
in our method. But given that the rate of two-body encoun- 
ters scales with the progenitor mass (see e.g. §7.5 of Binney 
& Tremaine 2008) both effects can be again minimized by 
applying the Minimum Entropy Method to tidal debris that 
originate from low-mass systems. 

Finally, it is well known that observational errors plus the 
presence of a background in the debris sample may also in- 
troduce biases in methods of statistical inference. For exam- 
ple, on account of our position within the Galaxy a systematic 
increase in the position and velocity errors with heliocentric 
distance seems unavoidable. The large number of ways in 
which errors and systematic biases may arise in phase-space 
catalogues prevents a generic approach to this problem, being 
more convenient to inspect this issue in individual cases. 

We plan to examine these and other aspects of our method 
in a forthcoming contribution, where realistic simulations 
of self-gravitating stellar clusters undergoing tidal disruption 
will be built by means of N-body techniques. Also, it will 
prove an exciting exercise to explore to what extent local 
gravity tests may complement those devised on cosmologi- 
cal scales (e.g. Zhao, Peacock & Li 2012) with help of mock 
Gaia catalogues. 
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